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Abstract 

Anisotropy  is  present  in  the  angular  distributions  of  neutrons  departing  from  a 
nuclear  scattering  event.  This  anisotropy  cannot  be  defined  in  a  closed- form  solution,  as 
in  the  Klein-Nishina  distribution  for  gamma  rays  following  scattering  events,  nor  is  the 
degree  and  behavior  of  anisotropy  only  dependent  on  the  incident  energy  of  the  particle. 
In  fact,  for  neutrons  leaving  a  scattering  event,  the  anisotropic  behavior  of  the  angular 
distribution  is  dependent  on  the  incident  energy  of  the  neutron,  the  type  of  scatter 
being  elastic  or  inelastic  along  with  the  inelastic  level,  and  the  species  struck. 

The  underlying  question  is,  if  anisotropic  behavior  is  worth  the  computational 
cost  to  be  included  in  certain  simulations,  and  if  so,  what  level  of  precision  is  effected  by 
the  inclusion  of  anisotropic  scatter. 

A  Watt  spectrum  of  U235  fission  neutrons  was  examined  as  it  collided  with  species 
in  a  nitrogen-oxygen  atmosphere.  In  a  stochastic  manner,  10'  collision  samples  were 
taken,  utilizing  cross  section-based  weighting  for  random  sampling  of  collision  types  and 
cross  section  weighting  along  with  concentration  weighting  to  determine  the  species 
struck.  The  collective  anisotropy  of  the  resultant  angular  distribution  was  apparent, 
with  a  definite  average  forward  bias  across  the  spectrum  and  a  bias  toward  scattering 
angles  less  than  30  degrees. 
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Additionally,  when  lower  energies  are  eliminated  from  consideration,  the  forward 
bias  increases.  This  leads  to  the  conclusion  that,  on  average,  the  higher  the  energy,  the 
greater  the  apparent  magnitude  of  the  forward  bias  of  the  anisotropic  scatter. 

Using  1-D  slab  geometry,  two  studies  were  conducted  exploring  the  relative  effect 
of  anisotropic  scatter  as  compared  to  isotropic  scatter  in  the  center  of  mass  reference 
frame.  The  maximum  relative  error  of  0.24%  was  observed  in  the  energy  dimension  and 
0.23%  in  the  time  dimension.  This  can  serve  as  a  first  approximation  for  more  complex 
problem  geometries  and  more  robust  scatter  mechanics.  In  short,  if  precision  is  required 
past  the  second  decimal  place  in  long-distance  high-altitude  transport  utilizing  isotropic 
scatter  in  the  center  of  mass  reference  frame,  anisotropy  in  the  center  of  mass  reference 
frame  deserves  consideration. 
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THE  EFFECT  OF  ANISOTROPIC  SCATTER  ON  ATMOSPHERIC 

NEUTRON  TRANSPORT 

I.  Introduction 

In  general,  the  response  of  any  object  or  medium  to  neutron  radiation  depends  on 
the  properties  of  the  object  and  energy  of  that  particle.  Any  system’s  response  also 
depends  on  the  angle  of  incidence  of  the  neutron.  Combining  the  energy  and  angle  of 
incidence,  in  the  rest  frame  of  the  object  or  system  of  concern,  dictates  the  response  the 
system  will  have  to  that  incident  neutron.  This  object  of  concern  might  be  a  sensor  of 
some  type  or  particles  in  the  medium  subjected  to  or  exposed  to  the  neutron  flux. 

The  response  of  a  system  can  change  based  on  the  quantity  of  neutrons,  not  just 
the  overall  energy  deposited  by  a  neutron  flux.  More  plainly,  a  small  number  of  high 
energy  neutrons  as  compared  to  a  large  number  of  low  energy  neutrons  will  generate 
different  responses,  even  though  the  sum  of  energy  might  be  the  same. 

Any  model  for  this  type  of  neutron  transport  will  incorporate  an  energy  range  of 
interest,  times  of  arrival  of  neutrons,  and  flux  or  fluence.  Additionally,  based  on  the 
problem  parameters  and  available  input  data,  a  given  model  will  have  requirements  for 
precision  and  accuracy.  Some  applications  may  tolerate  large  errors  and  others  may 
need  extremely  precise  calculations.  This  precision  and  accuracy  requirement  drives  the 
level  of  fidelity  required  within  a  model. 
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Application 


The  Air  Force  Technical  Applications  Center  (AFTAC)  monitors  compliance 
with  nuclear  test  ban  treaties.  Reliability  of  the  U.S.  Nuclear-detonation  Detection 
System  (USNDS)  requires  observing  multiple  phenomenologies,  including  neutron 
radiation.  AFIT  is  developing  a  new,  more-capable  Monte  Carlo  neutron  radiation 
transport  simulator  in  support  of  these  efforts. 

HASTEN 

High-Altitude  and  Space  Transport  Estimator  -  Neutrons  or  HASTE-N  is  a 
software  package  being  developed  by  Dr.  Kirk  Mathews.  HASTE-N,  in  general,  is 
designed  to  simulate  high  fidelity,  high  altitude  neutron  transport  using  stochastic 
Monte-Carlo  methods  in  an  effort  to  generate  accurate  modeling  for  specific  types  of 
low-density  medium  transport  problems.  [1] 

The  cornerstone  of  HASTE-N,  and  why  it  is  different  than  general  purpose 
stochastic  transport  codes,  is  based  on  the  inclusion  of  physical  processes  that  are 
normally  not  included  in  neutral  particle  transport.  One  key  process  is  the  independent 
motion  of  the  material  in  the  problem,  including  the  neutron  source,  the  transport 
medium,  and  the  object  of  interest  itself.  Typically  these  are  all  stationary,  which  is  not 
the  case  in  high  altitudes  with  detectors  in  orbit.  Another  key  process  is  the  influence 
of  gravity  on  the  neutron.  For  most  neutral  particle  transport  problems  in  short  time 
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scales  or  over  small  distances,  gravity  is  not  included.  Finally,  the  thermal  motion  of 


atoms  in  the  transport  medium  is  considered  in  order  to  more  accurately  predict  the 
energy  of  a  neutron  following  a  scatter.  [2] 

The  HASTE-N  code  is  also  developed  to  take  advantage  of  some  of  the  special 
features  of  the  problem  set  of  long  distance  low-density  medium  transport, 
predominantly  long  distances  between  collisions,  changing  medium  density,  and 
gravitational  effects.  Properties  of  low  density  air  as  the  medium  and  known 
gravitational  parameters  allow  optimization  of  the  code  for  tasks  like  cross-section 
lookup,  mean  free  path  calculations,  root-solving  problems. 

Additionally,  HASTE-N  allows  the  user  to  toggle  the  inclusion  of  different 
physical  processes  independently.  This  permits  individual  processes  to  be  examined  for 
verification  purposes  both  during  development  and  by  end  users.  An  individual  process 
can  then  be  studied  for  its  net  effect  on  a  given  problem.  The  added  benefit  of  this 
feature  is  the  ability  to  trade  fidelity  for  execution  speed  when  required. 

HASTE-N-TE 

HASTE-N  Test  Environment  (HASTE-N-TE)  is  a  software  package  developed  by 
Captain  Whitman  Daily  to  provide  a  simulation  environment  for  the  testing  of  HASTE- 
N  modules  separate  from  the  complexity  of  the  overall  HASTE-N  software  package. 
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For  this  reason,  HASTE-N-TE  functions  as  the  environment  used  to  test  the  modules 


discussed  in  this  document. 

Development  Focus 

The  underlying  question  is  if  anisotropic  behavior  is  worth  the  computational 
cost  to  be  included  in  certain  HASTE-N  problem  runs,  and  if  so,  what  level  of  precision 
is  effected  by  the  inclusion  of  anisotropic  scatter.  In  order  to  address  this  question,  first 
anisotropic  scatter  was  implemented,  then  verified,  and  finally,  the  measurement  of  the 
net  effect  on  an  overall  simulation. 

The  focus  of  this  document  is  on  the  development  of  a  small  part  of  the  HASTE- 
N  software  package  as  a  whole,  namely  the  modules  covering  scatter  mechanics  and 
cross-section  look  ups.  The  effort  described  here  was  to  implement  and  verify  these 
portions  of  HASTE-N  individually  to  contribute  to  the  verification  of  the  software 
package  as  a  whole. 

The  first  module  covers  the  anisotropic  scatter  of  neutrons  following  a  scattering 
interaction  with  a  species  of  nuclei  in  the  air.  In  general,  implementation  of  anisotropic 
behavior  consists  of  interpreting  and  interpolating  Legendre  polynomial  coefficients  from 
ENDF  files,  constructing  probability  density  functions  and  cumulative  distribution 
functions  from  interpolated  coefficients,  and  evaluating  the  angular  distribution  to 
determine  scattering  angle. 
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The  second  module  is  focused  on  the  ingestion  and  use  of  cross-section  data  from 


ENDF  file  extractions.  These  cross  sections  along  with  atmospheric  composition  data 
are  used  by  HASTE-N  to  determine  the  species  struck,  mean  free  path  length,  the  type 
of  interaction,  and  the  inelastic  level  of  scattering  events.  Additionally,  because  cross- 
sections  are  used  so  frequently  during  calculations  of  each  particle  history,  emphasis  was 
placed  on  execution  speed  of  cross-section  look  ups. 

Verification 

Following  the  implementation  and  testing  of  the  anisotropic  scatter  and  cross- 
section  modules,  a  stripped  down  version  of  HASTE-N  referred  to  as  HASTE-N  Test 
Environment  (HASTE-N-TE)  developed  by  Capt.  Whitman  Daily  was  used  to  verify 
and  test  the  behavior  of  the  two  modules.  Additionally,  parameter  studies  were 
conducted  to  observe  the  effect  of  anisotropy  on  long  distance  high  altitude  transport 
problems  as  simulated  by  HASTE-N-TE. 
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II.  Theory 


The  Anisotropic  Scatter  of  Neutrons 

Following  collision  with  a  nucleus,  an  elastic  or  inelastic  collision,  the  neutron 
departs  the  nucleus  with  a  given  deflection  angle  from  the  initial  direction  of  travel. 

This  angular  distribution  is  uniform  about  the  initial  direction  of  travel  of  the  neutron, 
but  the  deviation  from  that  initial  direction  is  drawn  from  a  probability  distribution 
function  and  varies  from  0  to  n  radians  (or  0  to  180  degrees). 

Consider  a  neutron,  initially  travelling  in  direction  fl'  =  v  j v' ,  that  collides  with 

a  nucleus  and  leaves  the  collision  in  a  new  direction  fl .  The  cosine  of  the  deflection 
angle  for  the  collision  is  usually  denoted  as  fiQ  =  iY  •  D .  This  is  the  cosine  of  a  polar 
angle  for  spherical  coordinates  with  the  pole  aligned  with  ST  .  The  other  spherical 
coordinate  angle,  representing  rotation  around  that  axis  is  conventionally  denoted  as  luq 
.  For  isotropic  scattering,  p()  is  be  uniformly  distributed  between  -1  (backscat.ter)  and 
+1  (forward  scatter),  and  is  uniformly  distributed  between  — 7T  and  7T .  The 
subscript  zero  is  dropped  for  conciseness  of  notation  in  the  rest  of  this  thesis. 

The  classical  mechanics  solution  for  the  scattering  collision  of  two  spheres  is  that 
the  distribution  in  isotropic  in  the  center  of  mass  (CM)  frame.  Thus,  the  probability 
density  function  (PDF)  is 


6 


[2.1] 


f isotropic  (^) 
^ “isotropic  (^) 


0  else 

V.  — 7T  <  U)  <  +7T 

/  Z7T 

0  else 


Using  this  classical  approximation,  if  the  mass  of  the  nucleus  is  much  greater  than  that 
of  a  neutron,  such  as  for  neutron  scattering  in  Uranium,  it  suffices  for  many  applications 
to  approximate  the  angular  distribution  as  isotropic  in  the  lab  frame.  However,  in 
scattering  a  neutron  from  an  atom  of  air,  the  mass  ratio  is  from  about  14:1  for  nitrogen- 
14  to  about  18:1  for  oxygen  18.  In  these  cases,  isotropic  in  the  lab  frame  is  not  an 
acceptable  approximation  to  isotropic  scatter  in  the  CM  frame.  The  relation  between 
the  scattering  cosine  in  the  lab  frame  and  in  the  CM  frame  for  elastic  scatter,  is  simply 


V CM  ~ 


[2.2] 


where  A  is  the  ratio  of  the  mass  of  the  nucleus  to  that  of  the  neutron.  For  elastic  scatter 
with  A  >  1 ,  the  PDF  in  the  lab  frame  is 


A  Cartesian  plot  of  the  PDF  for  isotropic  scatter,  for  A=14  and  A=18,  is  presented  in 
figure  1. 
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Figure  1:  Probability  density  function  for  isotropic  scatter  in  the  lab  frame,  plotted  in  Cartesian 
coordinates.  Curves  are  displayed  for  A— 14,  .<4=18,  and  A=co. 

The  steeper  curve  (blue)  is  the  PDF  for  the  less  massive  atom,  A  =  14 .  The  red 
curve  is  for  A  =  18  .  The  dashed  line  is  the  PDF  if  scattering  were  isotropic  in  the  lab 
frame,  i.e.,  for  yl  =  oo  .  It  is  easier  to  visualize  the  extent  of  the  anisotropy  using  a  plot 
of  the  surface  of  the  distribution  function  f  ^/j,  ^ ,  oj j .  Figure  2  shows  this  function  for  a 
neutron  that  was  initially  moving  directly  to  the  right.  Thus,  the  pole  on  the  right  of 


the  not-quite-sphere  is  forward  scatter  and  the  one  on  the  left  is  backscatter. 


Figure  2:  A  plot  of  the  surface  of  the  distribution  function  f  for  a  neutron  incident  from 

the  left  moving  right.  Note  the  plot  is  rotationally  symmetric  about  the  direction  of  travel  of  the 
incident  neutron. 

However,  these  plots  are  always  rotationally  symmetric  around  the  axis,  so  one 
usually  finds  plots  that  are  a  cross-section  of  this  shape  where  it  intersects  with  a  plane 
through  the  origin  perpendicular  to  the  line  of  sight.  Such  a  plot  is  a  polar  parametric 
plot  with  n  as  the  independent  variable,  /(h)  as  the  radius  and  angular  coordinate 
9  =  arctanl/i  .  Figure  3  presents  the  distributions  from  figure  2  in  this  way.  Forward 
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scatter  is  for  p  =  1 ,  whence  0  =  0,  which  is  to  the  right  in  figure  3.  Sideways  scatter  is 
for  p  =  0  ,  whence  9  =  tv  /  2  .  This  plot  makes  it  easy  to  see  the  extent  of  anisotropy. 
Remember  that  this  is  a  plot  of  the  PDF  of  p ,  not  the  PDF  of  9 . 


Figure  3:  A  polar  parametric  plot  with  p  as  the  independent  variable,  /  (p)  as  the  radius  and 
angular  coordinate  0  =  arctan(p).  Curves  are  displayed  for  A=14,  d=18,  and  A= oo. 

This  type  of  angular  distribution  is  commonly  encountered  in  the  content  of 

photons  deflecting  off  free  electrons,  resulting  in  the  Klein-Nishina  formula.  In  the  case 
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of  photons  colliding  with  electrons,  this  distribution  is  only  dependent  on  the  initial 
energy  and  is  well-defined  in  a  closed-form  equation. 

The  remainder  of  this  thesis  examines  the  extent  to  which  scattering  in  the  CM 
frame  is  not  actually  isotropic  (due  to  quantum-mechanical  effects).  In  all  subsequent 
formulas  and  plots,  p  denotes  the  scattering  deflection  cosine  in  the  center  of  mass 
frame. 


The  actual  anisotropic  angular  distribution  of  neutrons  (in  the  CM  frame) 
following  a  collision  with  a  nucleus,  is  dependent  on  the  initial  energy  of  the  collision  in 
the  center  of  mass,  the  species  of  nucleus  struck  (both  atomic  number  and  weight),  and 
the  excited  level  of  the  nucleus  as  a  result  of  an  inelastic  collision.  Unfortunately,  there 
is  no  known  closed  form  solution  for  this  angular  distribution.  Experimentally,  these 
distributions  are  measured.  ENDF  presents  coefficients  for  truncated  Legendre 
polynomial  series  that  are  used  to  approximate  the  distribution  functions  in  the  form 

j= 0 


Legendre  Polynomials 

Legendre  polynomials,  in  general,  are  azimuthally  symmetric  solutions  to 
Legendre’s  differential  equation, 
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+  n  (n  +  l)Pn  (a;)  =  0 


[2.5] 


A 

dx 


1-ad 


—  P 

dx  n 


where  P;  i^xj  is  the  nth  Legendre  polynomial.  Solving  the  above  ordinary  differential 
equation  is  not  related  to  the  approximation  of  an  angular  distribution  function,  but  the 
fact  that  the  polynomial  solutions  to  the  ODE  are  orthogonal  to  each  other  makes  these 
polynomials  ideal  for  the  estimation  of  a  function  that  does  not  have  a  closed  form;  the 
use  of  additional  polynomials  from  the  generated  series  allows  achievement  of  the 
desired  precision  of  the  estimator  function. 


The  Legendre  polynomials  from  n=  0  to  n= 5  are  as  follows: 


po 

(x) 

=  1 

pl\ 

=  X 

P2 

(x) 

=Xi 

&2-1) 

P3 

(x) 

=xi 

5x3  -  3ar) 

P4 

(x) 

-  V\ 
/  8  1 

(35x4  -  30x2  +  3 ) 

P5 

(x) 

-  Vi 
/  8  ' 

63a;5  -  70a;3  +  15a; 

V 

The  polynomials  from  equation  [2.6]  are  shown  graphically  in  Figure  4  below. 
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Figure  4:  Graph  of  the  First  Five  Legendre  Polynomials.  The  polynomials  are  shown  on 
the  interval  ranging  (-1,1)  and  are  shown  absent  any  weighting  coefficients. 

(htt  p :  // en.  wikipedia.  org  /  wiki /Legendre  _  polynomials) 

Legendre  Polynomial  Generation 

The  Legendre  polynomials  are  typically  defined  and  generated  recursively  using 
Bonnet’s  recursion  formula: 


PQ(x 

PAx 


P  \x 

n 


1 

X 

(2 n  -  l^P^Q)  -  (n  -  l)Pw_2(a;) 
n 


[2.7] 
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Though  able  to  be  generated  through  recursion  and  expanded  into  single 
polynomials  as  in  equation  [2.6],  the  polynomials  can  be  evaluated  during  the  recursion, 
resulting  in  less  possibility  of  computational  error  due  to  repeated  arithmetic  operations. 


Using  Legendre  Polynomials 

These  polynomials  are  used  here  as  an  estimator  of  f(p) ,  the  probability  density 
function  (PDF)  of  the  cosine  of  the  scattering  angle  of  a  neutron  and  nucleus  collision. 
The  polynomials,  individually  weighted  and  summed  together,  form  the  function 
estimator,  /(h)-  [3] 


n 

Y  (1]PM) 

3= 0 

0  else 


[2.8] 


where 


a .  = 
3 


2.7  + 1 


a . 

3 


[2.9] 


and  the  a  coefficients  are  tabulated  in  ENDF.  To  save  execution  time  in  the  Monte 

3 

Carlo  code,  the  coefficients  a.  are  computed  and  stored  in  arrays  as  the  coefficients  a . 
are  read  in  from  the  ENDF  files. 


The  factor  UUU- 
2 


arises  in  the  orthogonalization  of  the  Legendre  polynomials: 
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[2.10] 


£rijlpMpMiii=ei*  = 

A  PDF  is  not  only  normalized  to  integrate  to  1,  but  also  must  be  non- negative 
everywhere.  The  former  is  guaranteed  for  «(|  =  ;  the  latter  depends  on  the  values  of 

the  remaining  coefficients.  Presuming  non-negativity,  it  can  be  integrated  to  produce 
the  CDF, 


1  j  =  k 
0  j  v*  k 


f(m)  = ,  [2-11] 

which  must  be  monotone  non-decreasing. 

Cumulative  Distribution  Function  Evaluation 

In  order  to  retrieve  a  useful  value  from  the  PDF  in  equation  [2.8],  as  in  the 
chance  that  a  neutron  scatters  in  a  specific  direction  range,  the  PDF  must  be  integrated 
for  use  as  a  CDF.  The  CDF  is 


n-T  1 


fM=  FHx) dx  =  ^  i212i 


3  =  0 


j= 0 


where  the  coefficients  b .  are  given  by 


b.  = 


co  ci 

C  1  —  C  .  .  ! 
3~l  J+ 1 


3  =  0 

j  =  l,..,n  —  1 
j  =  n,  n  +  1 
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[2.13] 


with 


c . 


j 


2j  +  l 


[2.14] 


Like  the  PDF,  the  CDF  can  be  generated  and  evaluated  recursively.  This  is  done  for 
two  reasons,  the  first  for  the  minimizing  of  arithmetic  error  propagation  and  the  second 
for  an  increase  in  execution  speed. 


Neutron  Cross  Sections 

In  general,  a  cross  section  governs  the  probability  that  a  given  interaction  can 
occur.  In  the  context  here,  the  probability  that  an  incident  neutron  will  interact  with  a 
nucleus,  sometimes  referred  to  as  total  cross  section  for  interaction,  &totai-  Several 
interactions  of  a  neutron  with  a  nucleus  can  occur,  including  absorption,  elastic  scatter 
(n,n),  inelastic  scatter  (n,n’),  and  other  less  common  types. 

Cross  Section  for  Air 

Cross  section  data  is  readily  available  for  dry  air  at  sea  level  at  standard 
temperature  and  pressure  and  is  well-supported  by  experimental  data.  For  use  in 
HASTE-N,  this  data  would  not  accurately  depict  the  problem,  since  the  software  is 
designed  to  define  the  atmosphere’s  composition  for  each  problem.  A  more  consistent 
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approach  is  to  use  the  defined  atmosphere  to  construct  a  total  cross  section  for  the 
defined  atmosphere. 


For  example,  for  total  cross  section  for  interaction  for  air,  O’tota^air 


l 

^  _  V"'  Asotope,air 

u  total{ air)  Z  j  •’k  to£aZ(isotope  k) 

k=l 
n 

_ V"'  f  element, air  Asotope, element 

/  j  /  j  Jj  Jj  utoia/(element  i, isotope  j) 

t  If  1 


[2.15] 


where  l  is  the  total  number  of  species  in  the  atmosphere,  n  is  the  number  of  elements  in 
the  atmosphere,  m.  is  the  number  of  isotopes  of  element  i,  /“otope,a*r  is  the  atmospheric 
fraction  of  isotope  k,  and  so  on. 


Based  upon  how  the  atmospheric  concentrations  are  specified,  &totauair\  can  be 
constructed  in  two  ways,  both  displayed  in  equation  [2.15].  If  each  species  to  be 
considered  in  the  atmosphere  (elements  or  isotopes)  is  listed  along  with  the  fraction  of 
the  air  it  comprises,  the  cross  sections  of  each  species,  weighted  by  their  atmospheric 
fraction  (/))  are  summed  to  produce  the  total  cross  section.  If  the  atmosphere  is  instead 
further  specified  by  the  fraction  of  component  elements  (/) ,  along  with  the  fraction  of 
occurrence  of  isotopes  of  that  element  (/[) ,  the  cross  sections  of  each  species  is  weighted 
by  both  the  elemental  composition  of  the  atmosphere  along  with  the  natural  occurrence 
of  the  isotopes  of  that  element,  then  summed  to  produce  the  total  cross  section. 
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This  total  cross  section  for  interaction  for  air  can  then  be  used  for  mean  free  path 


calculations  within  HASTE-N 


The  Watt  Spectrum 

To  define  any  fission  spectrum  in  the  Evaluated  Nuclear  Data  Library,  ENDL  [4], 
the  Watt  spectrum  can  be  used, 

W  (a,  b,E^J  =  C  e~aE  sinh  (JbE^j  [2.16] 


where  E  is  the  energy  of  the  emitted  neutron  and  the  normalization  constant,  C,  is 
obtained  by  integration  of  W: 


C  = 


1 


sinh  ( \fbE )  dE 


[2.17] 


The  coefficients  a  and  b  are  different  for  each  individual  isotope  as  well  as  slight 
variances  based  on  incident  neutron  energy. 

To  approximate  the  neutron  fission  spectrum  for  U235  as  caused  by  neutrons  in 
the  MeV  range,  a  =  0.7  MeV-1  and  b  =  1.0  MeV-1 .  [5]  This  is  not  a  perfect 
approximation  for  the  produced  neutron  spectrum  of  U235,  but  will  be  used  as  a 
reasonable  baseline  for  the  purposes  of  parameter  studies  in  this  document.  Further 
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details  on  use  of  the  Watt  spectra  for  fission  spectra  can  be  found  in  UCRL-TR-203351. 
[4]  (Note:  That  report  contains  an  incorrect  formula  for  the  normalization  constant.) 
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III.  Methodology 


Since  HASTE-N  is  being  developed  in  Fortran,  the  same  development 
environment  was  selected  for  implementation  of  the  scatter  angle  distribution  and  cross 
section  modules.  The  functions  and  subroutines  were  written  to  interface  with  the 
existing  development  effort,  and  utilized  existing  routines  from  other  HASTE-N 
modules.  Microsoft  Visual  Studio  2013  along  with  the  Intel  Math  Kernel  Library  was 
used  to  develop,  debug,  and  test  the  modules. 

Angular  Distribution  Implementation 

Since  the  departure  angle  of  a  neutron  following  a  scatter  event  is  dependent  on 
not  just  the  species  of  nuclei  struck  and  the  incident  energy  of  the  neutron  in  the  center 
of  mass  frame,  but  also  on  the  type  of  scatter  (elastic  or  inelastic),  a  large  data  set  of 
Legendre  coefficients  is  required  for  each  isotope.  For  each  isotope,  the  ENDF  database 
has  varying  energy  levels  for  elastic  collisions  as  well  as  varying  energy  levels  for  the 
first  thirty  inelastic  levels.  One  of  the  design  parameters  of  this  module  was  not  to  use 
a  proprietary  data  file  format;  files  directly  from  ENDF  can  be  read  by  the  software 
without  alteration,  allowing  HASTE-N  to  be  updated  simply  by  adding  newer  ENDF 
files  as  they  are  released. 
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Legendre  Coefficient  Data 


Dynamically  sized,  multidimensional  data  structures  were  constructed  to  hold  the 
Legendre  coefficients.  Conceptually,  the  data  is  organized  in  a  four  dimensional  matrix 
by  isotope,  energy,  inelastic  level,  and  the  Legendre  coefficient  subscript.  At  the  time  of 
ingestion  of  the  ENDF  files,  every  one  of  those  four  dimensions  must  be  determined 
dynamically. 

The  number  of  isotopes  is  based  on  the  complexity  of  the  atmosphere  for  a  given 
problem.  Elemental  composition  of  the  atmosphere  as  well  as  isotopes  of  each  element 
are  considered,  based  on  the  level  of  detail  required  for  the  atmosphere’s  composition  in 
a  specific  problem.  As  a  reasonable  first  approximation,  the  initial  atmospheric 
composition  was  selected  to  contain  the  five  most  common  species,  in  order  of 
commonality:  N14,  N15,  01(>,  O18,  and  O1'.  The  fraction  of  each  element  can  be  specified 
as  well  as  the  fraction  of  each  isomer  of  each  element  in  problem  parameters. 

In  a  given  ENDF  file  containing  Legendre  Coefficients,  coefficients  for  multiple 
incident  energies  are  listed.  The  number  of  energies  listed  is  primarily  dependent  on 
how  much  data  was  gathered  experimentally.  [6]  For  elastic  collisions,  the  listed  incident 
energies  are  from  as  low  as  ~1.0  x  10  eV  and  can  range  into  the  GeV  level  depending 
on  the  measurements  recorded.  For  the  inelastic  level  data,  the  lowest  energy  listed  is 
the  excitation  energy  for  that  inelastic  excitation  level,  or  the  Q.  For  the  purposes  of 
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this  software,  the  incident  energy  data  extends  no  higher  than  14.5  MeV,  but  this  can 


be  raised,  provided  the  provided  ENDF  data  extends  to  the  desired  energy  level  bound. 

ENDF  provides  data  for  elastic  collisions  and  for  the  first  thirty  inelastic  levels  of 
most  isotopes,  then  a  continuum  for  the  levels  above  thirty.  Fortunately,  with  the 
upper  bound  for  energy  set  at  14.5  MeV,  only  the  inelastic  levels  with  Q  values  less  than 
the  upper  energy  bound  need  be  considered.  The  software  dynamically  determines  the 
levels  required  for  each  isotope  based  on  the  specified  maximum  incident  energy. 

The  four-dimensional  Legendre  data  structure  contains  the  Legendre  coefficients 
themselves.  Based  on  the  experimental  data  collected  and  the  degree  of  anisotropy 
present  for  a  given  isotope,  energy,  and  inelastic  level  combination,  anywhere  from  2  to 
32  coefficients  may  be  provided  for  energies  below  the  upper  energy  bound  of  14.5  MeV. 
[6]  More  than  32  coefficients  may  be  provided  for  higher  energies,  and  the  software  can 
easily  be  modified  for  this  application  if  desired.  It  should  be  noted  that  larger 
magnitudes  and  greater  numbers  of  coefficients  indicate  a  more  intense  anisotropy  of  the 
selected  scatter  and  hence  a  more  complex  function  describing  that  distribution. 

Once  the  collision  type  (isotope,  energy,  and  inelastic  level)  has  been  fully 
determined  through  weighted  random  sampling,  a  linear  interpolation  of  the  coefficients 
from  the  energy  above  and  below  the  target  energy  produces  the  coefficients  used  for  the 
PDF  and  the  CDF. 
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Cross  Section  Implementation 


Once  a  collision  location  has  been  determined  using  mean  free  path  calculations, 
cross  sections  of  all  the  possible  collision  types  are  used  to  determine  the  species  struck 
and  the  collision  type.  During  this  weighted  random  sampling,  certain  collision  types 
can  be  suppressed,  such  as  absorption  or  less  common  interactions  like  (n,2n)  reactions. 
To  determine  the  collision  type,  two  weighted  random  samples  must  be  executed;  the 
first  sampling  selects  the  species  struck  and  the  second  sampling  determines  the  type  of 
scattering  event  (elastic  or  inelastic)  and  the  inelastic  level. 

To  determine  the  species  struck,  total  cross  sections  for  each  species  in  the 
atmosphere,  weighted  by  the  atmospheric  fraction  of  each  species,  need  to  be  considered. 
Additionally,  in  order  to  utilize  embedded  random  number  generators  in  the  Math 
Kernel  Library  without  additional  mathematical  operations,  the  random  sample  has 
been  constructed  to  utilize  random  numbers  on  their  intrinsic  range,  R  G  [0, 1) . 

The  test  used  to  select  the  struck  species  is 

k 

E£  g.total 
■'li 

R  <  — -  [3.11 

n 

Er  ptotal 

Hu  i 
i= 1 
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where  R  is  a  pseudo-random  number,  is  the  total  cross  section  for  interaction  for 
the  'i,u  species  in  the  atmosphere  and  f  is  the  fraction  of  the  atmosphere  represented  by 
the  species,  and  n  is  the  total  number  of  species  in  the  atmosphere.  The  test  begins 
with  k=l  and  is  repeated  with  the  same  random  number,  incrementing  k  until  the  test 
returns  true  or  k=n- 1.  If  the  test  is  true,  the  kth  isotope  is  struck,  and  if  k=n- 1  the  nth 
isotope  is  known  to  be  struck  without  performing  the  final  test.  To  minimize  the 
number  of  arithmetic  operations,  in  successive  tests,  results  from  the  previous 
summations  can  be  used,  as  well  as  ordering  the  list  of  species  with  the  highest 
atmospheric  concentrations  first 

Once  the  species  struck  has  been  determined,  the  type  of  scatter  (elastic  or 
inelastic  with  level)  must  be  determined.  First,  the  maximum  inelastic  level  is  selected 
by  comparing  the  energy  of  the  incident  neutron  in  the  inertial  frame  of  the  struck 
nucleus  to  the  excitation  energies  (Q)  of  the  available  inelastic  levels.  Any  level  with  a 
Q  greater  than  the  incident  energy  (in  the  frame  of  the  center  of  mass  of  the  collision) 
cannot  be  reached  and  is  therefor  excluded.  In  a  similar  manner  to  equation  [3.1],  the 
level  selection  test  is 


^ scatter 

R  <  1^0 - 

n 

^ scatter 

i= 0 


[3.2] 
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where  ascatte 7  is  the  cross  section  for  scatter  and  excitation  to  the  d 1 1  inelastic  level  with 

l 

i=  0  being  elastic,  and  n  is  the  total  number  of  available  inelastic  levels  based  on  the  Q 
value.  The  test  begins  with  k=  0  (elastic  scatter)  and  is  repeated  with  the  same  random 
number,  incrementing  k  until  the  test  returns  true  or  k=n- 1.  If  the  test  is  true,  the 
scatter  is  the  kth  inelastic  level,  and  if  k=n- 1  the  scatter  takes  place  at  the  nth  inelastic 
level,  the  highest  level  possible  for  that  incident  neutron  energy  without  performing  the 
final  test.  The  summations  occur  during  each  iteration  of  the  test,  but  summation 
results  from  previous  iterations  can  be  used  to  minimize  arithmetic  operations.  In  all 
examined  cases,  for  all  energies  under  14.5  MeV,  the  elastic  cross  section  is  greatest,  and 
hence  the  most  likely  scatter  is  elastic.  In  that  elastic  case,  only  one  test  occurs.  Even 
though  this  test  is  O (n),  functionally  it  is  much  faster  due  to  the  frequency  of  elastic 
scatter. 

Unified  energy  grid 

Since  so  many  cross  section  lookups  are  required  to  fully  define  one  particle’s 
history,  including  mean  free  path  calculations  in  air  and  multiple  collision  calculations 
as  in  equations  [3.1]  and  [3.2],  a  unified  energy  grid  is  created  prior  to  running  a 
problem  set  in  an  effort  to  speed  up  individual  cross  section  look  ups. 

In  a  given  ENDF  cross  section  file,  the  energies  are  not  standardized  across  the 
whole  ENDF  library,  but  are  based  on  round  numbers  (like  1.0e5  eV  or  3.5e-6  eV)  and 
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on  excitation  energy  levels  for  the  particular  species  being  measured.  Consequently, 
each  ENDF  cross  section  file  has  cross  sections  for  a  different  set  of  energies. 

The  unified  energy  grid  is  a  three  dimensional  matrix;  the  first  dimension  is 
energy,  the  second  is  the  cross  section  type,  and  the  third  is  the  species  or  isotope. 

Much  like  the  Legendre  coefficient  data,  the  size  of  two  of  these  three  dimensions  must- 
be  determined  dynamically,  based  on  the  problem  definition  and  the  contents  of  the 
ENDF  libraries  being  used. 

To  construct  the  size  of  the  first  dimension,  every  energy  listed  in  every  ENDF 
file  for  the  defined  atmosphere  is  read  in  and  placed  in  a  list.  Duplicates  are  eliminated 
from  the  list,  and  the  list  is  sorted  and  truncated  based  on  upper  and  lower  energy 
bounds.  This  incorporates  all  the  energies  of  interest  for  all  of  the  isotopes  taken  into 
consideration. 

The  cross  section  types  included  as  the  second  dimension  are  the  cross  section  for 
elastic  scatter,  the  cross  section  for  total  inelastic  scatter,  and  the  total  cross  section  for 
interaction.  Finally  the  third  dimension  is  comprised  of  all  the  species  included  in  the 
defined  atmosphere. 

Once  the  dimensions  of  the  grid  are  determined,  the  cross  sections  themselves 
need  to  be  populated  into  the  unified  grid.  Some  energy  /  cross  section  pairs  directly 


map  to  the  ENDF  files,  but  some  pairs  in  the  grid  will  not  directly  map  to  the  provided 
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data  and  need  to  be  interpolated.  Linear  interpolation  is  used  to  avoid  introducing 
arithmetic  precision  errors,  as  multiple  linear  interpolations  are  mathematically 
equivalent  to  a  single  linear  interpolation.  This  is  not  the  case  with  some  other 
interpolation  methods. 

Once  the  grid  is  constructed,  look  ups  into  the  grid  will  frequently  not  land  on  an 
energy  boundary  and  will  require  interpolation,  also  linear.  So,  the  cross  section  value 
used  for  interpolation  may  be  interpolated  twice,  once  when  the  unified  grid  is 
populated,  and  once  when  a  table  look  up  occurs.  Since  linear  interpolating  twice  to  a 
value  is  equivalent  to  linear  interpolating  once  to  the  same  value,  no  interpolation  error 
is  introduced. 


Linear  Interpolation 

Given  a  standard  linear  interpolation  of  a  function  f(x)  with  known  values  at  x0 

and  xt, 


f(x)  =  f(x0)  +  [f(x1)-f(xQ 


Xl~Xo) 


If  another  interpolated  point,  x\  is  calculated  in  the  same  interval, 


f(xo)  +  {f(xl)-f(x  0 


X, 


xn 


[3.3] 


[3.4] 
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it  can  be  shown  that  x  can  be  found  by  interpolating  x'  with  either  x{)  or  xv 


f(x)  =  f(xQ)  +  (f(x 


I  \ 


[3.5] 


By  substituting  equation  [3.4]  into  equation  [3.5], 


f(x) 
f  W 
f(x) 


=  flx  n  + 


fM  +  ifM-flx, 


=  f{x0)  +  (f(x1)-f(x0) 
=  f(xo)  +  (f(xi)-f(xo) 


X1  xo. 


X1  xo 


-f(x  I 


[3.6] 


equation  [3.3]  is  produced,  showing  multiple  linear  interpolations  are  mathematically 
equivalent  to  a  single  linear  interpolation. 


Verification  Design 

As  with  most  program  development  efforts,  the  developed  modules  are  tested 
independently  from  the  larger  programming  effort,  in  order  to  focus  the  search  for 
compilation,  run  time,  and  logic  errors  in  the  code.  The  majority  of  the  verification 
work  was  performed  in  Microsoft  Visual  Studio  2013  [7]  and  graphical  verification  was 
performed  in  Matlab  2014b  [8]. 
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Anisotropy  Testing 


The  first  question  to  be  addressed  with  regards  to  anisotropic  behavior  of 
neutrons  following  scattering  events  is,  to  what  degree  is  anisotropic  behavior  present 
even  in  the  simplest  of  conditions?  With  the  provided  ENDF  data,  prior  to  any 
simulation  effort,  the  anisotropic  angular  distribution  for  the  considered  species  was 
examined  to  determine  if  the  anisotropic  behavior  was  evident  at  all. 

Collision  Testing 

With  the  use  of  the  Watt  spectrum  for  neutrons  produced  from  the  fission  of  U235 
along  with  the  atmospheric  concentrations  specified  for  a  five  species  atmosphere  (N14, 
N15,  O16,  O18,  and  O17),  the  collision  choosing  functionality  was  tested.  This  serves  not 
only  to  ensure  all  possible  collision  types  were  represented  and  functioning  correctly,  but 
also  serves  as  an  indicator  of  the  general  likelihood  of  individual  scatter  events  occurring 
during  a  more  expansive  simulation. 

Net  Effect  of  Anisotropy 

Finally,  as  an  examination  of  the  effect  of  anisotropy  of  neutron  scatter  on 
problems  of  this  type,  two  parameter  studies  were  conducted  using  HASTE-N-TE, 
incorporating  atmospheric  species  concentrations,  weighted  scattering  occurrence  rates, 
and  energies  generated  by  a  line  source.  Both  of  these  studies  utilized  a  flat-earth, 
stationary  atmosphere  model  varying  altitude  and  energy  respectively. 
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IV.  Results  and  Analysis 


A  Mat-lab  2014b  script  was  written  to  construct  the  PDF  from  Legendre 
weighting  coefficients  combined  with  their  associated  polynomials  and  to  plot  them  as  in 
figure  3.  Because  the  incident  particles  here  have  mass,  the  energy  is  in  the  center  of 
mass  frame  of  the  collision,  and  may  differ  from  the  energy  of  the  incident  neutron  in 
the  lab  reference  frame  based  on  the  movement  of  the  nucleus  being  struck. 

Anisotropic  Angular  Distributions 

Prior  to  examination  of  any  results  from  simulations,  the  first  step  was  to 
examine  the  anisotropic  behavior  of  neutrons  in  general,  using  coefficients  directly  from 
ENDF,  to  determine  any  patterns  or  common  behavior  of  the  angular  distributions  in 
general.  Note:  The  areas  under  different  PDFs  plotted  as  in  figure  1  are  all  equal  to  1, 
by  the  normalization  of  a  PDF.  However,  the  areas  inside  of  different  PDFs  on 
parametric  polar  plots,  as  in  figure  3  and  in  this  chapter,  are  not  related  to 
normalization  and  can  be  very  different. 

First,  the  behavior  of  the  angular  distribution  of  a  neutron  in  the  exit  channel 
following  an  elastic  collision  was  considered.  The  isotope  was  fixed,  and  only  energy 
was  varied.  The  first  isotope  examined  in  this  manner  was  O16,  and  energy  levels  of  1 
MeV,  4  MeV,  and  7  MeV  are  shown  in  figure  5  below. 
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Figure  5:  /(mcm)  plotted  radially  against  0(nCM)  of  Neutrons  at  1  MeV,  4  MeV,  and  7 
MeV  following  elastic  scattering  events  with  O16.  These  energies  are  measured  in  the 
reference  frame  of  the  center  of  mass  of  the  collision.  With  changing  energies  of  incident 
neutron  colliding  with  the  same  species  in  only  elastic  scattering  events,  significantly 
different  angular  distributions  result. 

Knowing  that  isotropic  behavior  would  be  displayed  as  a  circle  of  radius  0.5, 
anisotropic  scatter  is  clearly  present,  as  seen  by  two  distinct  lobes  in  the  1  MeV 
distribution  and  three  distinct  lobes  in  the  4  MeV  and  7  MeV  distributions.  The 
changing  number  of  lobes  shows  the  behavior  is  not  just  scaled  based  on  energy,  but  can 
have  distinct  features  appear  and  disappear  along  a  given  energy  range.  Another 
important  feature  to  note  is  that  at  1  MeV,  there  is  no  forward  scatter  or  backscatter 
bias,  but  this  is  not  true  at  higher  energies. 
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The  same  plot  in  Figure  5  was  repeated  in  Figure  6,  only  changing  the  species 
struck  to  N14.  Only  elastic  scatter  is  represented,  and  1  MeV,  4  MeV,  and  7  MeV  are 
represented. 


Figure  6:  /(mcm)  plotted  radially  against  0(nCM)  of  Neutrons  at  1  MeV,  4  MeV,  and  7 
MeV  following  elastic  scattering  events  with  N14.  These  energies  are  measured  in  the 
reference  frame  of  the  center  of  mass  of  the  collision.  With  changing  energies  of  incident 
neutron  colliding  with  the  same  species  in  only  elastic  scattering  events,  significantly 
different  angular  distributions  result. 

Similar  to  Figure  5,  the  1  MeV  distribution  shows  no  forward  or  backscatter  bias, 
but  the  forward  scatter  bias  is  clearly  evident  at  higher  energies.  Again,  the  lobes  of  the 
distributions  appear  to  present  themselves  and  disappear  as  energy  changes. 
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Once  it  was  clear  that  anisotropic  scatter  is  present  across  energy  levels,  inelastic 


level  scatter  was  examined.  Since  it  is  now  established  that  the  anisotropic  behavior 
can  vary  with  energy,  the  next  task  would  be  to  determine  if  anisotropic  behavior  varies 
with  the  level  of  an  inelastic  scatter. 
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Figure  7:  /  plotted  radially  against  0(nCM)  of  8  MeV  Neutrons  following  2nd,  4th,  and 
6th  level  inelastic  scattering  events  with  O18.  These  energies  are  measured  in  the  reference 
frame  of  the  center  of  mass  of  the  collision.  With  the  same  energy  of  incident  neutron 
colliding  with  the  same  species,  based  on  the  inelastic  level,  significantly  different  angular 
distributions  result. 

Figure  7  shows  the  angular  distribution  for  varying  levels  of  anisotropy  across 
inelastic  levels,  where  the  angular  distributions  from  three  inelastic  level  scatters  (2nd, 
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4th,  and  6th)  from  8  MeV  neutrons  are  incident  to  0lb.  Three  distinct  scattering 
behaviors  are  observed.  First,  in  the  2nd  level  inelastic  scatter  four  lobes  are  present, 
indicating  a  forward  scatter,  backscatter,  and  direct  side  scatter  bias  with  zeros  between 
those  lobes.  Next,  in  the  4th  level  inelastic  scatter  a  slight  forward  scatter  and 
backscatter  bias  is  observed,  with  smooth  transitions  throughout.  Finally,  in  the  6th 
level  inelastic  scatter,  a  slight  bias  to  side  scatter  is  observed,  again  with  smooth 
transitions. 

Figures  5,  6  and  7  are  not  present  to  show  patterns,  or  to  draw  definitive 
conclusions,  but  they  do  clearly  show  that  anisotropic  behavior  is  present,  and  that  the 
behavior  can  vary  across  any  of  the  parameters  of  the  collision,  namely  species  struck, 
inelastic  level,  and  incident  neutron  energy.  However,  there  is  no  expectation  that  a 
clear  pattern  would  emerge.  If  there  was  a  distinct  pattern,  a  closed-form  solution 
would  likely  be  common  practice  for  determining  neutron  scatter  behavior. 

Instead,  the  conclusion  that  can  be  drawn  is  that  anisotropic  behavior  is  clearly 
present  in  neutron  scatter  events,  and  bears  investigation,  as  it  may  have  significant 
effect  on  certain  simulations  related  to  long-distance,  point  source  problems. 

Weighted  Random  Choosing 

Next,  in  order  to  verify  both  the  cross  section  method  for  choosing  collision  types 
and  the  Watt  spectrum  for  determining  incident  neutron  energies,  100  million  energies 
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were  drawn  from  the  Watt  spectrum,  and  collisions  types  were  tabulated  using  cross 


sections  and  atmospheric  concentrations. 

During  this  process  of  weighted  random  sampling,  individual  neutrons  with  their 
associated  energies  were  examined,  and  virtually  every  possible  collision  type  and  species 
was  encountered,  ensuring  that  the  chance  of  each  type  of  collision  was  appropriately 
weighted  as  shown  in  Equation  [3.1]  and  in  Equation  [3.2], 
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Table  1:  Listing  of  Collisions  Produced  by  108  Neutrons  Drawn  from  the  Watt  Spectrum. 
Atmospheric  species  fractions  along  with  cross  sections  were  used  in  weighting  the 
collisions. 


Level 

N14 

N15 

o16 

o18 

o17 

0 

72886617 

387202 

25451238 

7492 

5388 

1 

213075 

1142 

13427 

1529 

629 

2 

439105 

426 

89258 

186 

41 

3 

76752 

298 

18133 

132 

77 

4 

165195 

178 

8093 

331 

27 

5 

43304 

117 

2775 

144 

11 

6 

74840 

128 

576 

69 

34 

7 

24413 

21 

484 

76 

12 

8 

11900 

20 

597 

25 

17 

9 

31891 

5 

307 

64 

0 

10 

8992 

14 

354 

59 

9 

11 

5236 

21 

295 

20 

3 

12 

9230 

6 

198 

33 

1 

13 

2761 

7 

243 

31 

2 

14 

251 

5 

9 

3 

4 

15 

181 

7 

4 

8 

2 

16 

227 

4 

10 

2 

2 

17 

1869 

2 

2 

2 

2 

18 

2817 

5 

2 

2 

3 

19 

2865 

4 

0 

5 

1 

20 

1836 

0 

0 

6 

2 

21 

728 

0 

0 

3 

1 

22 

1182 

0 

0 

7 

6 

23 

898 

0 

0 

5 

0 

24 

658 

0 

0 

2 

0 

25 

470 

0 

0 

1 

1 

26 

274 

0 

0 

2 

0 

27 

180 

0 

0 

3 

0 

28 

114 

0 

0 

1 

2 

29 

0 

0 

0 

2 

0 

30 

0 

0 

0 

0 

0 
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Table  1  functions  as  a  display  of  the  collision  choosing  portion  of  the  module. 


Several  elements  of  the  table  are  present  as  expected,  such  as  the  majority  of  collisions 
are  incident  to  N14,  the  most  plentiful  species  in  the  defined  atmosphere.  Additionally, 
the  majority  of  the  collisions  are  elastic  (98.77%),  represented  by  the  0th  level. 

The  Effect  of  Anisotropy 

With  the  individual  elements  of  the  module  functioning  in  accordance  with 
theory,  the  underlying  question  is  whether  anisotropic  behavior  is  worth  the 
computational  cost  to  be  included  in  certain  HASTE-N  problem  runs,  and  if  so,  what 
level  of  precision  is  effected  by  the  inclusion  of  anisotropic  scatter.  The  first  part  of  this 
question  is  addressed  here,  and  the  latter  will  become  apparent  as  the  HASTE-N 
software  package  goes  through  further  testing  and  development. 

First,  given  the  reasonable  set  of  parameters  for  a  sample  problem  listed  in  Table 
2  (coefficients  for  the  Watt  spectrum  representing  neutrons  emitted  from  the  fission  of 
U235,  and  a  five  species  atmosphere  consisting  of  two  isotopes  of  nitrogen  and  three 
isotopes  of  oxygen  with  associated  concentrations),  10b  collisions  were  selected  by 
random  sampling  producing  the  collisions  listed  in  Table  1. 
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Table  2:  Problem  Parameters 


Atmospheric  Fractions  Watt  Spectrum  Coefficients 


Nitrogen 

0.78084 

a 

0.7 

Oxygen 

0.20946 

b 

1 

Isotopic  Concentrations 


N14 

0.99636 

N15 

0.00364 

016 

0.99757 

018 

0.00205 

017 

0.00038 

Next,  the  energy  of  each  incident  neutron  causing  each  of  the  10b  events  was  used 
to  determine  the  species  struck  and  the  inelastic  level.  This  information  was  then  used 
to  interpolate  tabulated  Legendre  coefficients.  These  10s  sets  of  coefficients  were 
averaged,  producing  an  average  PDF  of  first  scatters,  f  ,  for  the  problem  as  a  whole, 
representing  the  overall  scatter  angle  distribution  for  the  problem  parameters,  taking 
into  account  weighting  from  collision  type,  incident  energy,  and  species  struck.  This 
overall  f  ^fxj  is  graphed  in  Figure  8  below. 
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Figure  8:  /(Mcm)  plotted  radially  against  0(/j,cm)  of  10s  neutrons  drawn  from  the  Watt 
spectrum  for  U235  fission  incident  to  an  atmosphere  as  specified  in  Table  2.  The 
anisotropic  nature  of  the  angular  distribution  PDF  is  apparent,  with  a  significant  forward 
bias. 

The  overall  f  y/ij  does  not  have  any  distinct  lobes  or  sharp  features,  as  it  is 
averaged  over  10b  different  function  estimators.  Readily  apparent  is  a  forward  bias, 
with  a  distinct  bias  for  deflection  angles  less  than  30  degrees.  This  shows  that  the 
problem  set  as  a  whole  has  significant  anisotropic  behavior  and  should  be  considered  in 
problems  with  parameters  similar  to  those  in  Table  2. 
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Because  lower  energies  are  shown  to  show  less  forward  bias,  as  seen  in  Figure  5 


and  Figure  6,  the  same  averaging  was  examined  excluding  lower  energies.  Considering 
this  anisotropic  behavior,  excluding  lower  energies  could  be  useful  if  a  given  detection 
system  has  a  lower  energy  threshold,  or  if  due  to  the  geometry  of  the  problem,  which 
has  not  been  considered  here,  lower  energy  neutrons  need  not  be  considered. 


Figure  9:  f  (hCm)  plotted  radially  against  0(fj,CM)  in  varying  energy  ranges  of  108  neutrons 
drawn  from  the  Watt  spectrum  for  U235  fission  incident  to  an  atmosphere  as  specified  in 
Table  2.  Distributions  shown  suppress  lower  energies.  When  higher  energies  are 
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displayed  absent  lower  energies  the  foreword  bias  is  more  apparent,  showing  the  forward 

bias  comes  from  the  higher  energy  neutrons  more  than  the  lower  energy  neutrons. 

When  the  same  10s  neutrons  are  considered,  masking  out  lower  energies,  the 
foreword  bias  becomes  more  apparent.  In  fact,  as  a  larger  range  of  low  energies  is 
masked,  as  in  Figure  9,  bias  is  more  prominent. 

This  leads  to  the  conclusion  that,  on  average,  the  higher  the  energy,  the  greater 
the  apparent  magnitude  of  the  forward  bias  of  the  anisotropic  scatter.  Based  on  the 
design  of  the  problem,  including  geometry  and  detector  sensitivity,  this  further  increased 
anisotropy  of  scatter  angle  could  have  an  even  greater  effect  on  the  fidelity  of  the  model. 

Altitude  Parameter  Study 

Through  the  use  of  HASTE-N-TE,  a  study  varying  altitude  while  examining 
anisotropy  was  conducted.  HASTE-N-TE  was  configured  to  utilize  a  14.06  MeV 
isotropic  line  source,  forcing  a  leakage-suppressed  first  elastic  scatter  only,  and  then 
populating  energy-time  bins  through  the  use  of  a  point  estimator. 

The  atmosphere  used  in  the  study  is  of  uniform  density  and  is  composed  of 
isotopic  concentrations  listed  in  Table  2  and  stretches  from  0  to  86  km  above  the  earth 
with  vacuum  above.  The  atmosphere  is  1-D  slab  geometry,  above  a  flat  earth.  The 
detector  was  placed  at  30,000  km  and  the  source  location  was  varied  in  5  km  increments 
from  5  km  to  85  km.  Using  these  parameters,  two  runs  of  108  particles  at  each  altitude 
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were  simulated,  the  first  run  using  isotropic  scatter,  then  the  second  using  anisotropic 


scatter  as  determined  by  the  scattering  module. 

Prior  to  examining  any  results,  knowing  the  line  source  is  14.06  MeV  and  that 
nitrogen- 14  is  the  most  common  species  in  the  atmosphere,  the  anisotropic  behavior  of 
14  MeV  neutrons  elastically  scattered  off  nitrogen- 14  was  examined  and  is  displayed  in 
Figures  10  and  11  below. 


120 
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Figure  10:  /  (hcm)  plotted  radially  against  d(/j,CM)  of  14  MeV  neutrons  following  elastic 
scattering  events  with  N14  including  detail  of  side  and  backscatter  regions.  This 
distribution  is  measured  in  the  reference  frame  of  the  center  of  mass  of  the  collision  and 
generated  directly  from  ENDF  coefficients.  Scaled  to  arbitrary  units  of  intensity. 
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Figure  11:  /  (/zCM)  plotted  against  [iCM  of  14  MeV  neutrons  following  elastic  scattering 
events  with  N14,  units  of  |i  vs.  arbitrary  intensity  units  of  the  generated  PDF,  f(p).  This 
distribution  is  measured  in  the  reference  frame  of  the  center  of  mass  of  the  collision  and 
generated  directly  from  ENDF  coefficients. 

Examining  Figures  10  and  11,  a  forward  positive  bias  presents  itself  at  scatter 
angles  less  than  27°  with  a  negative  bias  present  for  side  and  back  scatter  regions.  The 
uniform  intensity  of  0.5  representing  isotropic,  the  point  where  the  forward  bias  end  can 
be  determined  by  setting  f(p)  equal  to  0.5. 


/  0)  =  0.5 

p  =  cos(27°)  =  0.893 


[4.1] 
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Even  though  the  lobes  representing  side  scatter  shown  in  Figure  10  are 


significantly  less  than  isotropic,  it  should  be  noted  that  the  center  of  a  lobe  in  the  polar 
plot  of  the  distribution  can  be  mathematically  located  by  setting  the  derivative  of  f(p) 
equal  to  zero. 

/'O)  =  o 

H  =  cos(85°)  =  0.087  ' 

In  order  to  examine  the  time-energy  bins  from  each  of  the  runs  of  this  study,  to 
best  measure  the  relative  effect  of  anisotropy,  the  time  and  energy  dimensions  were 
examined  separately.  First,  neutron  counts  in  time-integrated  energy  bins  were  plotted 
at  varying  altitudes,  then  neutron  counts  in  energy-integrated  time  bins  were  plotted  at 
the  same  varying  altitudes.  Collapsing  these  dimensions  allows  a  good  display  of  the 
relative  difference  of  anisotropic  scatter  as  compared  to  isotropic  scatter  in  the  time  and 
energy  dimensions. 


44 


Relative  Error 
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Figure  12:  Energy  bin  neutron  counts  vs.  relative  error  over  five  source  altitudes. 

Relative  error  is  calculated  with  anisotropic  treated  as  the  true  value  and  isotropic 
treated  as  the  measured  value. 

The  relative  error  between  neutron  counts  in  isotropic  and  anisotropic  time- 
integrated  energy  bins, 


isotropic 


-V, 


anisotropic 


' rel 


anisotropic 


[4.3] 
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as  shown  in  Figure  12,  has  several  distinct  features  that  are  directly  linked  to 
anisotropic  scatter  behavior.  Thus,  where  the  relative  error  is  positive,  a  calculation 
using  isotropic  scatter  in  the  CM  frame  overestimates  the  result  of  a  calculation  that 
uses  the  ENDF  anisotropic  scatter  in  the  CM  frame. 


First,  in  a  direct  backscatt.er  collision  event,  where  the  neutron  deposits  the 
maximum  amount  of  energy  in  the  atom,  and  has  minimal  exit  channel  energy,  the 
isotropic  scatter  model  overestimates  counts  by  nearly  0.24%.  This  occurs  where  the 
difference  between  the  isotropic  distribution  and  the  anisotropic  distribution  of  nitrogen- 
14  diverge  the  most,  at  a  direct  backscatt.er. 


The  maximum  energy  loss  in  a  collision,  A =  E'  —  E  ,  for  a  14  MeV  neutron, 
E'  =  14  MeV ,  is 


A  E 


max 

loss 


=  E' 


4  m  m 

_ n  a 

(m  +  m  ) 

\  n  a) 


«  3.5  MeV  , 


[4.4] 


so  that  the  minimum  energy  of  the  neutron  after  the  backscat.ter  //  =  —  1 ,  is 


E  .  =  14  MeV  —  A E  pa  10.5  MeV . 

mm  max 

loss 


[4.5] 
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This  is  seen  in  Figure  12.  It  should  be  noted  that  in  equation  [4.4]  the  change  in 
reference  frame  is  not  accounted  for,  as  exit  channel  energy  remains  in  the  center  of 
mass  frame  of  the  collision  in  this  study. 

Second,  an  analogous  and  opposite  behavior  presents  itself  in  the  region  of 
forward  scatter,  where  energy  transfer  is  minimal,  and  exit  channel  energy  is  near  14 
MeV.  Here,  an  underestimate  by  isotropic  scatter  is  evident,  a  result  of  the  distinct- 
forward  bias  of  nitrogen- 14  anisotropic  scatter  in  the  forward  region,  as  seen  in  Figure 
10. 


The  dip  in  the  center  of  the  graph  in  Figure  12  is  a  combination  of  both  the 
design  parameters  of  the  simulation,  and  of  the  sideways  scatter  bias  of  oxygen-16,  the 
second  most  common  species  in  the  defined  atmosphere.  Here,  due  to  the  1-D  slab 
geometry  and  flat-  earth  model  used,  some  neutrons  travel  nearly  parallel  t-o  the  top  of 
the  atmosphere  for  a  long  distance  before  scattering,  then,  the  point  estimator  only 
attenuates  the  neutron  slightly  before  it  exits  the  atmosphere  and  is  tallied.  This  is 
more  pronounced  at-  higher  source  elevations,  as  the  source  is  closer  t-o  the  top  of  the 
defined  atmosphere.  Oxygen- 16  exhibits  an  elastic  side  scatter  bias  at  higher  energies, 
which  is  why  the  isotropic  scatter  model  is  an  underestimate. 

At  the  highest  source  elevation,  85  km,  some  noise  is  present  at  the  higher  energy 
bins,  or  the  bins  populated  by  forward  scatter  events  because  the  scatter  is  forced  to 
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occur  in  a  small  amount  of  the  atmosphere,  and  when  assigned  the  correspondingly  low 


weight  using  Monte-Carlo  methods,  produces  poor  counting  statistics. 


Energy-Integrated  Time  Bins 


Figure  13:  Time  bins  neutron  counts  vs.  relative  error  over  five  source  altitudes.  Relative 
error  is  calculated  with  anisotropic  treated  as  the  true  value  and  isotropic  treated  as  the 
measured  value. 

Now,  examining  the  same  time-energy  bin  relative  error  data  as  in  Figure  12,  but 
with  energy-integrated  energy  bins  shown  in  Figure  13,  several  additional  features  can 
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be  noted.  First,  the  relative  error  is  zero  in  figure  13,  between  0.625  and  0.64  seconds 
depending  on  source  altitude.  This  corresponds  to  a  scattering  angle  of  27°.  At  this 
angle,  the  anisotropic  angular  distribution  for  nitrogen-14  intersects  with  the  isotropic 
distribution  case,  as  seen  in  Figure  10.  This  behavior  is  expected,  at  the  angle  where  an 
anisotropic  scatter  distribution  matches  the  value  of  an  isotropic  distribution  the  error 
between  the  two  would  be  zero,  and  the  expected  difference  of  the  time-energy  bins 
populated  by  that  particular  scatter  should  also  be  near  zero. 

Second,  at  -0.66  seconds,  a  significant  peak  in  the  error  is  present,  getting  larger 
at  higher  source  altitudes.  This  phenomena  is  a  direct  result  of  the  negative  bias 
toward  backscatter  of  nitrogen- 14.  The  isotropic  scatter  model  significantly 
overestimates  the  weight  of  a  particle  going  down  into  the  atmosphere  from  the  source 
then  backscatt.ering  to  the  detector.  It  follows  that,  for  a  higher  source  location,  there  is 
more  atmosphere  to  travel  through,  and  hence  the  collision  is  given  a  higher  weight, 
leading  to  increased  relative  error  with  increasing  source  altitudes.  This  leads  to  a 
maximum  of  0.23%  relative  error. 

Finally,  at  times  past  0.7  seconds,  a  significant  increase  in  noise  is  present  and  is 
not  consistent  with  increasing  altitudes.  This  is  primarily  an  artifact  of  the  1-D  slab 
geometry  combined  with  the  backscatter  mechanics  of  nitrogen-14  and  oxygen-16. 

One  way  to  achieve  a  flight  time  of  0.75  seconds  is  as  follows, 
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[4.6] 


fc4MeV>(0.75-fW,  I^^MeV) 

v  m  ^  ’  1  y  m 

V  n  V  n 

d L2  +  d22  =~  300002 
t  =  0.71  sec 
d1  =  1446km 
d2  —  29965km 

where  dl  is  the  distance  traveled  prior  to  the  scatter,  and  d2  is  the  distance  between  the 
scatter  and  the  detector.  Here,  the  neutron  travels  nearly  parallel  to  the  flat  earth  for 
1446  km  and  then  up  to  the  detector.  This  is  only  possible  with  a  1-D  slab  geometry 
flat  earth  model.  Since  scattering  events  similar  to  the  example  have  an  extremely  low 
weight  due  to  the  long  distance  traveled  within  the  atmosphere,  they  produce  poor 
counting  statistics,  hence  the  noise. 

The  noise,  however,  is  not  consistently  positive  or  negative  relative  error.  This  is 
as  a  result  of  the  species  struck.  Nitrogen-14,  the  most  common  species  struck  in  the 
atmosphere,  has  a  negative  side-scatter  bias  and  oxygen-16,  the  second  most  common 
species,  has  a  positive  side-scatter  bias.  This,  combined  with  poor  counting  statistics 
means  single  events  have  a  significant  effect  on  the  noise  depending  on  the  species 
struck,  nitrogen-14  showing  isotropic  producing  an  overestimate,  and  oxygen-16 
producing  an  underestimate  for  the  scattering  events  that  populate  time  bins  in  this 
region. 
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Energy  Parameter  Study 


Through  the  use  of  HASTE-N-TE,  a  study  varying  energy  while  examining 
anisotropy  was  conducted.  HASTE-N-TE  was  configured  to  utilize  11  isotropic  line 
source  energies  ranging  from  1  MeV  to  14  MeV,  forcing  a  leakage-suppressed  first  elastic 
scatter  only,  and  then  populating  energy-time  bins  through  the  use  of  a  point  estimator. 

The  atmosphere  used  in  the  study  is  of  uniform  density  and  is  composed  of 
isotopic  concentrations  listed  in  table  2  and  stretches  from  0  to  86  km  above  the  earth 
with  vacuum  above.  The  atmosphere  is  1-D  slab  geometry,  above  a  flat  earth.  This 
configuration  is  identical  to  that  used  in  the  altitude  parameter  study. 

The  detector  was  placed  at  30,000  km  and  the  source  location  was  fixed  at  43 
km,  the  center  of  the  defined  atmosphere.  Using  these  parameters,  two  runs  of  10s 
particles  at  each  energy  were  simulated,  the  first  run  using  isotropic  scatter,  then  the 
second  using  anisotropic  scatter  as  determined  by  the  scattering  module. 

Since  flight  times  vary  both  with  the  scatter  mechanics  and  the  energy  of  the  line 
source,  the  results  of  this  study  are  best  viewed  in  the  relative  error  of  time-integrated 
energy  bins.  This  allows  analysis  of  the  change  in  the  shape  of  the  relative  error  vs. 
energy  plot  shown  in  Figure  14  as  incident  energy  varies. 
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Figure  14:  Energy  vs.  relative  error,  energy  varying  from  4  MeV  to  14  MeV.  Relative 
error  is  calculated  with  anisotropic  treated  as  the  correct  value  and  isotropic  treated  as 
the  approximate  value. 

In  the  energy  range  depicted  in  Figure  14,  the  shape  of  each  relative  error  vs. 
energy  plot  is  similar  to  that  seen  in  Figure  12,  with  the  same  interpretation  for  the 
features  in  each  plot  as  described  in  the  altitude  study.  Primarily,  this  is  evidence  that 
the  behavior  of  the  relative  error  of  time-integrated  energy  bins  is  consistent  across  an 
energy  range,  not  just  at  14.06  MeV  as  analyzed  in  the  altitude  study. 


Additionally,  across  the  energy  range  of  4  MeV  to  14  MeV  as  displayed  in  Figure 


14,  the  relative  error  increases  linearly  with  energy.  This  behavior  of  the  relative  error 
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can  be  used  as  a  determining  factor  when  deciding  to  apply  anisotropy  to  an  energy 


range,  based  on  the  noise  level  of  a  particular  simulation  or  the  precision  desired. 

At  energies  beneath  4  MeV,  however,  the  shape  of  the  relative  error  plot  begins 
to  change  as  seen  in  Figure  f5.  The  maximal  error  shifts  from  the  left  of  each  plot, 
meaning  the  region  of  energy  bins  populated  by  backscatters,  to  the  right  side  of  the 
plots,  the  region  of  energy  bins  populated  by  the  forward  scatters. 


0.6  0.8  1  1.2  1.4  1.6  1.8  2  2.2 

Energy.  MeV 


Figure  15:  Energy  vs.  relative  error,  energy  varying  from  1  MeV  to  2  MeV.  Relative 
error  is  calculated  with  anisotropic  treated  as  the  correct  value  and  isotropic  treated  as 
the  approximate  value. 

Due  to  the  decreasing  relative  error  in  the  1  MeV  to  2  MeV  energy  range 
combined  with  the  effect  of  anisotropy  being  reduced  for  this  energy  range  as  seen  in 
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Figure  6,  the  artifacts  of  the  problem  parameters  may  eclipse  any  conclusions  that  could 


be  drawn  from  this  behavior.  Higher  fidelity  simulation  in  HASTE-N  focused  on  this 
energy  range  could  allow  more  definitive  conclusions  regarding  the  behavior  of  the 
relative  error  in  this  energy  region. 


Figure  16:  Probability  density  functions  for  elastic  scattering  of  nitrogen-14,  1  MeV  to  2 
MeV.  Plot  is  of  |x  vs.  normalized  intensity  of  each  PDF.  Legendre  coefficients  are 
extracted  from  ENDF  without  interpolation. 

The  behavior  of  the  angular  distribution  probability  density  functions  for 
nitrogen-14,  as  shown  in  Figure  16,  shows  oscillations  in  the  behavior  of  the  PDFs  that 
are  not  easily  explained.  This  does  not  clearly  explain  the  behavior  of  the  relative  error 
of  energy  bins  in  this  region,  leading  to  the  supposition  that  the  shift  in  the  shape  of  the 
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relative  error  plots  in  Figure  15  is  not  related  to  anisotropic  behavior  changes,  but  from 
a  different  aspect  of  the  simulation. 
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V. 


Conclusion 


Anisotropy  is  present  in  the  angular  distributions  of  neutrons  departing  from  a 
nuclear  scattering  event.  This  anisotropy  cannot  be  defined  in  a  closed-form  solution  as 
in  the  Klien-Nishina  distribution  for  gamma  rays  following  scattering  events,  nor  is  the 
degree  and  behavior  of  anisotropy  only  dependent  on  the  incident  energy  of  the  particle. 
In  fact,  for  neutrons  leaving  a  scattering  event,  the  anisotropic  behavior  of  the  angular 
distribution  is  dependent  on  the  incident  energy  of  the  neutron,  the  type  of  scatter 
being  elastic  or  inelastic  along  with  the  inelastic  level,  and  the  species  struck. 
Additionally,  this  dependence  cannot  be  predicted  using  a  closed  form  solution,  and 
polynomial  functions  have  been  developed  and  are  tabulated  in  ENDF  based  on 
experimental  data. 

The  function  used  to  describe  the  estimator  of  the  PDF  of  the  angular 
distribution  for  a  given  set  of  conditions  is  constructed  of  weighted  Legendre 
polynomials.  In  general,  the  larger  the  magnitude  of  the  weighting  coefficients,  the 
more  severe  the  anisotropy.  Unfortunately  the  shape  of  the  angular  distribution  is  not 
predictable  based  on  any  consistent  factor,  and  experimental  data  must  be  used  to 
construct  angular  distributions  for  use  in  simulation  models. 

Amidst  all  of  the  other  computational  demanding  factors  of  a  Monte  Carlo 
neutron  transport  code  including  problem  geometry,  fidelity  demands,  cross  section 
lookups,  and  others,  inclusion  of  anisotropic  scattering  behavior  requires  examination  to 
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determine  if  the  increase  in  computational  cost  adds  appreciably  to  the  fidelity  of  the 


answer,  and  to  what  degree  fidelity  is  increased. 

A  spectrum  of  U235  fission  neutrons,  as  generated  by  the  Watt  spectrum,  was 
examined  as  it  collided  with  species  in  a  nitrogen-oxygen  atmosphere  with  normal 
natural  concentrations  of  nitrogen  and  oxygen  isotopes.  In  a  stochastic  manner,  10b 
collision  samples  were  taken,  utilizing  cross  section-based  weighting  for  random 
sampling  of  collision  types  and  cross  section  weighting  along  with  concentration 
weighting  to  determine  the  species  struck.  The  collective  anisotropy  of  the  resultant 
angular  distribution  was  apparent,  with  a  definite  average  forward  bias  across  the 
spectrum  and  a  bias  toward  scattering  angles  less  than  30  degrees. 

Additionally,  when  lower  energies  are  eliminated  from  consideration,  the  forward 
bias  increases.  This  leads  to  the  conclusion  that,  on  average,  the  higher  the  energy,  the 
greater  the  apparent  magnitude  of  the  forward  bias  of  the  anisotropic  scatter.  This  is 
worth  considering  when  undertaking  simulation  design  utilizing  detectors  with  specified 
sensitivity  ranges. 

Using  1-D  slab  geometry,  two  studies  were  conducted  exploring  the  relative  effect 
of  anisotropic  scatter  as  compared  to  isotropic  scatter  in  the  center  of  mass  reference 
frame.  The  maximum  relative  error  of  0.24%  was  observed  in  the  energy  dimension  and 
0.23%  in  the  time  dimension.  This  can  serve  as  a  first  approximation  for  more  complex 
problem  geometries  and  more  robust  scatter  mechanics.  In  short,  if  precision  is  required 
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past  the  second  decimal  place  in  long-distance  high-altitude  transport  utilizing  isotropic 


scatter  in  the  center  of  mass  reference  frame,  anisotropy  in  the  center  of  mass  reference 
frame  deserves  consideration. 

Further  Work 

Additional  examination  regarding  the  specific  nature  of  the  anisotropy  is 
suggested,  particularly  within  a  simulation.  Once  HASTE-N  is  able  to  function  with 
and  without  implementation  on  anisotropic  scatter,  detailed  parameter  studies  should  be 
conducted  to  determine  the  gain  in  fidelity  realized  by  the  inclusion  of  anisotropic 
neutron  scatter. 
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